!/**********************************************************
! * Author        : lvxj
! * Email         : jing3704@126.com
! * Last modified : 2024-08-08 15:51
! * Description   : 
! * *******************************************************/
module fangzhu
   implicit none 
   include 'mpif.h'
        integer myid,nps
        integer ierr
        integer V_INLET,WALL,INTERIOR
        real(8)::velo_iter_err
        integer max_iter_count_v
        integer ::Xmax,Ymax,Zmax  !Xmax,Ymax..grid
        integer ::Yst,Yed,Zst,Zed
        integer npy,npz,idy,idz
        real(8) relaxu,relaxv,relaxw,relaxp
        real(8) :: h
        real(8) :: nu
        real(8) :: Ls
        real(8) :: dt,Uw,end_time
        integer::istep,steps_st
        integer::MaxStep
        integer:: NUMAX
        real(8) ::dx,dy,dz,d
        integer::x1,x2,a_in,b_in
        real(8)::delt
        integer::y1,y2,z1,z2
        real(8) ::Vp,Sf1,Sf2,Sf3 
        real(8),allocatable,target ::u(:,:,:),v(:,:,:),w(:,:,:),p(:,:,:)
        real(8),allocatable,target ::u_ave(:,:,:),v_ave(:,:,:),w_ave(:,:,:),p_ave(:,:,:)
        integer,allocatable,target ::flagd11(:,:,:),flagd12(:,:,:),flagd13(:,:,:)
        real(8),allocatable,target ::un(:,:,:),vn(:,:,:),wn(:,:,:),pn(:,:,:)
        real(8),allocatable,target ::pc(:,:,:),pci(:,:,:)
        real(8),allocatable ::ui(:,:,:),vi(:,:,:),wi(:,:,:),pi(:,:,:)

        real(8),allocatable ::uf(:,:,:,:),vf(:,:,:,:),wf(:,:,:,:),pf(:,:,:,:)
        real(8),allocatable,target ::mf(:,:,:,:)
        real(8),allocatable ::uf0(:,:,:,:),vf0(:,:,:,:),wf0(:,:,:,:),pf0(:,:,:,:)
        real(8),allocatable ::mf0(:,:,:,:)
        real(8),allocatable,target ::pcf(:,:,:,:)

        real(8),allocatable ::gradux(:,:,:),graduy(:,:,:),graduz(:,:,:)
        real(8),allocatable ::gradvx(:,:,:),gradvy(:,:,:),gradvz(:,:,:)
        real(8),allocatable ::gradwx(:,:,:),gradwy(:,:,:),gradwz(:,:,:)
        real(8),allocatable ::gradpx(:,:,:),gradpy(:,:,:),gradpz(:,:,:)
        real(8),allocatable ::gradpcx(:,:,:),gradpcy(:,:,:),gradpcz(:,:,:)

        real(8),allocatable ::gradufx(:,:,:,:),gradufy(:,:,:,:),gradufz(:,:,:,:)
        real(8),allocatable ::gradvfx(:,:,:,:),gradvfy(:,:,:,:),gradvfz(:,:,:,:)
        real(8),allocatable ::gradwfx(:,:,:,:),gradwfy(:,:,:,:),gradwfz(:,:,:,:)
        real(8),allocatable ::gradpfx(:,:,:,:),gradpfy(:,:,:,:),gradpfz(:,:,:,:)
        real(8),allocatable ::gradpcfx(:,:,:,:),gradpcfy(:,:,:,:),gradpcfz(:,:,:,:)

        real(8),allocatable,target ::a_u(:,:,:)
        real(8),allocatable,target ::b_u(:,:,:)
        real(8),allocatable,target ::a_un(:,:,:)
        real(8),allocatable,target ::a_us(:,:,:)
        real(8),allocatable,target ::a_uw(:,:,:)
        real(8),allocatable,target ::a_ue(:,:,:)
        real(8),allocatable,target ::a_ut(:,:,:)
        real(8),allocatable,target ::a_ud(:,:,:)

        real(8),allocatable,target ::a_v(:,:,:)
        real(8),allocatable,target ::b_v(:,:,:)
        real(8),allocatable,target ::a_vn(:,:,:)
        real(8),allocatable,target ::a_vs(:,:,:)
        real(8),allocatable,target ::a_vw(:,:,:)
        real(8),allocatable,target ::a_ve(:,:,:)
        real(8),allocatable,target ::a_vt(:,:,:)
        real(8),allocatable,target ::a_vd(:,:,:)

        real(8),allocatable,target ::a_w(:,:,:)
        real(8),allocatable,target ::b_w(:,:,:)
        real(8),allocatable,target ::a_wn(:,:,:)
        real(8),allocatable,target ::a_ws(:,:,:)
        real(8),allocatable,target ::a_ww(:,:,:)
        real(8),allocatable,target ::a_we(:,:,:)
        real(8),allocatable,target ::a_wt(:,:,:)
        real(8),allocatable,target ::a_wd(:,:,:)

        real(8),allocatable ::a_p(:,:,:)
        real(8),allocatable ::b_p(:,:,:)

        real(8),parameter :: standard=1.0e-6
        integer,parameter :: max_iter=10
        integer:: counter

        integer gridsnum!,nonzero
end module

!------------------------------------------------------------------
SUBROUTINE TIMER(T)
IMPLICIT NONE 
INCLUDE 'mpif.h'
DOUBLE PRECISION,INTENT(INOUT) :: T
INTEGER*4 :: VAL(8)
CALL DATE_AND_TIME(VALUES=VAL)
T=DBLE(VAL(8))*1.D-3+DBLE(VAL(7))+ DBLE(VAL(6))*60.D0+DBLE(VAL(5))*3600.D0
RETURN
END

subroutine init()
   use fangzhu
   implicit none

   integer i
   integer ncpus
   integer lny,lnz,resy,resz

   V_INLET=1
   WALL=3
   INTERIOR=4
   Uw=1.0

   velo_iter_err=1.0e-7
   max_iter_count_v=10

   relaxu=0.3
   relaxv=0.3
   relaxw=0.3
   relaxp=0.3

   MaxStep=50001
   steps_st=100000!240000
   Xmax=100!240
   Ymax=100!160
   Zmax=100!100

   NUMAX=100!MaxStep
   dt=0.001
   end_time=100!0.101
   nu=1.0/1000!0.0008 !0.001
   Ls=1!10.0 !0.1
   dx=1.0/Xmax!0.1!1.0/diameter!0.1!Ls/Xmax
   dy=dx
   dz=dx
   d=dx
   h=dx
   Sf1=dy*dz
   Sf2=dx*dz
   Sf3=dx*dy
   Vp=dx*dy*dz


      open(12,file='np.dat')
      read(12,*)ncpus
      read(12,*)npy
      read(12,*)npz
      close(12)
      if(ncpus.ne.nps)then
        print*,'Wrong procesors! please check!'
        stop
      endif
      if((npz*npy).ne.nps)then
        print*,'Wrong npz and npy! please check!'
        stop
      endif
      if(nps.ne.1)then
        if(myid==0) then
          print*,'erro! please set nps == 1 !!!!'
        endif
        call MPI_Barrier(MPI_COMM_WORLD,ierr)
        stop
      endif
      idy=myid/npz
      idz=myid-idy*npz
      if(mod(Ymax,npy)==0)then
        lny=Ymax/npy
        Yst=idy*lny+1
        Yed=(idy+1)*lny
      else
        lny=Ymax/npy
        resy=mod(Ymax,npy)
        if(idy<resy)then
          Yst=idy*(lny+1)+1
          Yed=(idy+1)*(lny+1)
        else
          Yst=idy*lny+resy+1
          Yed=(idy+1)*lny+resy
        endif
      endif
      if(mod(Zmax,npz)==0)then
        lnz=Zmax/npz
        Zst=idz*lnz+1
        Zed=(idz+1)*lnz
      else
        lnz=Zmax/npz
        resz=mod(Zmax,npz)
        if(idz<resz)then
          Zst=idz*(lnz+1)+1
          Zed=(idz+1)*(lnz+1)
        else
          Zst=idz*lnz+resz+1
          Zed=(idz+1)*lnz+resz
        endif
      endif
    gridsnum=Xmax*(Yed-Yst+1)*(Zed-Zst+1)
!      print*,'myid,idy,idz,Yst,Yed,Zst,Zed:',myid,idy,idz,Yst,Yed,Zst,Zed

end subroutine init

subroutine alloc_mem()
   use fangzhu
   implicit none

   integer i,j,k,l

   allocate(u(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),v(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),w(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
    &   p(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(u_ave(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),v_ave(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
    &   w_ave(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),p_ave(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(flagd11(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),flagd12(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
    &   flagd13(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(un(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),vn(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
    &   wn(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),pn(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(a_u(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),a_v(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),a_w(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(b_u(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),b_v(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),b_w(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(pc(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),pci(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(ui(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),vi(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
   &    wi(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),pi(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))

   allocate(uf(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),vf(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    wf(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),pf(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(mf(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(uf0(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),vf0(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    wf0(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),pf0(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(mf0(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(pcf(Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))

   allocate(gradux(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),graduy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
   &    graduz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(gradvx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),gradvy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
   &    gradvz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(gradwx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),gradwy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
   &    gradwz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(gradpx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),gradpy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
   &    gradpz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))
   allocate(gradpcx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),gradpcy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2),&
   &    gradpcz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2))

   allocate(gradufx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),gradufy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    gradufz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(gradvfx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),gradvfy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    gradvfz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(gradwfx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),gradwfy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    gradwfz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(gradpfx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),gradpfy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    gradpfz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))
   allocate(gradpcfx(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),gradpcfy(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3),&
   &    gradpcfz(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2,3))

   allocate(a_ue(Xmax,Yst:Yed,Zst:Zed),a_uw(Xmax,Yst:Yed,Zst:Zed),a_un(Xmax,Yst:Yed,Zst:Zed),&
   &    a_us(Xmax,Yst:Yed,Zst:Zed),a_ut(Xmax,Yst:Yed,Zst:Zed),a_ud(Xmax,Yst:Yed,Zst:Zed))
   allocate(a_ve(Xmax,Yst:Yed,Zst:Zed),a_vw(Xmax,Yst:Yed,Zst:Zed),a_vn(Xmax,Yst:Yed,Zst:Zed),&
   &    a_vs(Xmax,Yst:Yed,Zst:Zed),a_vt(Xmax,Yst:Yed,Zst:Zed),a_vd(Xmax,Yst:Yed,Zst:Zed))
   allocate(a_we(Xmax,Yst:Yed,Zst:Zed),a_ww(Xmax,Yst:Yed,Zst:Zed),a_wn(Xmax,Yst:Yed,Zst:Zed),&
   &    a_ws(Xmax,Yst:Yed,Zst:Zed),a_wt(Xmax,Yst:Yed,Zst:Zed),a_wd(Xmax,Yst:Yed,Zst:Zed))
   allocate(a_p(Xmax,Yst:Yed,Zst:Zed),b_p(Xmax,Yst:Yed,Zst:Zed))
  

  u_ave=0.0
  v_ave=0.0
  w_ave=0.0
  p_ave=0.0

  u=0.0
  v=0.0
  w=0.0
  p=0.0
  un=0.0
  vn=0.0
  wn=0.0
  pn=0.0
  a_u=0.0
  a_v=0.0
  a_w=0.0
  pc=0.0
  pci=0.0
  ui=0.0
  vi=0.0
  wi=0.0
  pi=0.0

  uf=0.0
  vf=0.0
  wf=0.0
  pf=0.0
  mf=0.0
  pcf=0.0


  uf0=0.0
  vf0=0.0
  wf0=0.0
  mf0=0.0
  pf0=0.0

  gradux=0.0
  graduy=0.0
  graduz=0.0
  gradvx=0.0
  gradvy=0.0
  gradvz=0.0
  gradwx=0.0
  gradwy=0.0
  gradwz=0.0
  gradpx=0.0
  gradpy=0.0
  gradpz=0.0
  gradpcx=0.0
  gradpcy=0.0
  gradpcz=0.0

  gradufx=0.0
  gradufy=0.0
  gradufz=0.0
  gradvfx=0.0
  gradvfy=0.0
  gradvfz=0.0
  gradwfx=0.0
  gradwfy=0.0
  gradwfz=0.0

  a_ue=0.0
  a_uw=0.0
  a_un=0.0
  a_us=0.0
  a_ut=0.0
  a_ud=0.0

  a_ve=0.0
  a_vw=0.0
  a_vn=0.0
  a_vs=0.0
  a_vt=0.0
  a_vd=0.0

  a_we=0.0
  a_ww=0.0
  a_wn=0.0
  a_ws=0.0
  a_wt=0.0
  a_wd=0.0

  a_p=0.0
  b_p=0.0

  flagd11=0
  flagd12=0
  flagd13=0




end subroutine alloc_mem


subroutine dealloc_mem()
   use fangzhu
   implicit none

   deallocate(u,v,w,p)
   deallocate(u_ave,v_ave,w_ave,p_ave)
   deallocate(flagd11)
   deallocate(flagd12)
   deallocate(flagd13)
   deallocate(un,vn,wn,pn)
   deallocate(a_u,a_v,a_w)
   deallocate(b_u,b_v,b_w)
   deallocate(pc,pci)
   deallocate(ui,vi,wi,pi)

   deallocate(uf,vf,wf,pf)
   deallocate(mf)
   deallocate(mf0)
   deallocate(pcf)

   deallocate(gradux,graduy,graduz)
   deallocate(gradvx,gradvy,gradvz)
   deallocate(gradwx,gradwy,gradwz)
   deallocate(gradpx,gradpy,gradpz)
   deallocate(gradpcx,gradpcy,gradpcz)

   deallocate(gradufx,gradufy,gradufz)
   deallocate(gradvfx,gradvfy,gradvfz)
   deallocate(gradwfx,gradwfy,gradwfz)
   deallocate(gradpfx,gradpfy,gradpfz)
   deallocate(gradpcfx,gradpcfy,gradpcfz)

   deallocate(a_ue,a_uw,a_un,a_us,a_ut,a_ud)
   deallocate(a_ve,a_vw,a_vn,a_vs,a_vt,a_vd)
   deallocate(a_we,a_ww,a_wn,a_ws,a_wt,a_wd)
   deallocate(a_p,b_p)



end subroutine dealloc_mem

subroutine set_bdy()
   use fangzhu
   implicit none

   integer i,j,k,l
  l=3
  if(Zed==Zmax)then
     k=Zmax+1
     do j=Yst-1,Yed+1
    if(j>0.and.j<=Ymax)then
     do i=1,Xmax
        vf(i,j,k,l)=0.0!Uw !v_far
        mf(i,j,k,l)=0.0
        uf(i,j,k,l)=Uw !0.0 !u_far
        wf(i,j,k,l)=0.0 !w_far
     enddo
    endif
     enddo
  endif

  do k=Zst-1,Zed+1
    if(k>0.and.k<=Zmax)then
  do j=Yst-1,Yed+1
    if(j>0.and.j<=Ymax)then
  do i=1,Xmax
     flagd11(i,j,k)=INTERIOR
  enddo
    endif
  enddo
    endif
  enddo

  if(Zed==Zmax)then
     k=Zmax+1
     do j=Yst-1,Yed+1
    if(j>0.and.j<=Ymax)then
     do i=1,Xmax
        flagd11(i,j,k)=V_INLET
     enddo
    endif
     enddo
  endif

  if(Zst==1)then
     k=0
     do j=Yst-1,Yed+1
     do i=1,Xmax
        flagd11(i,j,k)=WALL
     enddo
     enddo
  endif
 
  i=0
  do k=Zst-1,Zed+1
  do j=Yst-1,Yed+1
     flagd11(i,j,k)=WALL
  enddo
  enddo

  i=Xmax+1
  do k=Zst-1,Zed+1
  do j=Yst-1,Yed+1
     flagd11(i,j,k)=WALL
  enddo
  enddo

  if(Yst==1)then
     j=0
     do k=Zst-1,Zed+1
     do i=0,Xmax+1
        flagd11(i,j,k)=WALL
     enddo
     enddo
  endif

  if(Yed==Ymax)then
     j=Ymax+1
     do k=Zst-1,Zed+1
     do i=0,Xmax+1
        flagd11(i,j,k)=WALL
     enddo
     enddo
  endif

  flagd12=flagd11
  flagd13=flagd11


end subroutine set_bdy

subroutine set_init_flux()
   use fangzhu
   implicit none

   integer i,j,k,l

   integer i1,i2,j1,j2

end subroutine set_init_flux

subroutine update(var)
   use fangzhu
   implicit none
   real(8)::var(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)

   integer south,north,top,down
   integer src,dst

   real(8)::buffer_sdy(0:Xmax+1,Zst:Zed,2),buffer_rvy(0:Xmax+1,Zst:Zed,2)
   integer i,j,k
   integer,dimension(MPI_STATUS_SIZE)::stat

   south=idy-1
   north=idy+1
   top=idz+1
   down=idz-1

   do j=Yst,Yst+1
   do k=Zst,Zed
   do i=0,Xmax+1
      buffer_sdy(i,k,j-Yst+1)=var(i,j,k)
   enddo
   enddo
   enddo

   if(south>=0)then
        dst=idz+south*npz
   else
        dst=MPI_PROC_NULL
   endif
   if(north<npy)then
        src=idz+north*npz
   else
        src=MPI_PROC_NULL
   endif

   call MPI_Sendrecv(buffer_sdy(0,Zst,1),(Xmax+2)*(Zed-Zst+1)*2,MPI_DOUBLE_PRECISION,dst,9000,&
   &    buffer_rvy(0,Zst,1),(Xmax+2)*(Zed-Zst+1)*2,MPI_DOUBLE_PRECISION,src,9000,MPI_COMM_WORLD,stat,ierr)


  if(north<npy)then
   do j=Yed+1,Yed+2
   do k=Zst,Zed
   do i=0,Xmax+1
      var(i,j,k)=buffer_rvy(i,k,j-Yed)
   enddo
   enddo
   enddo
  endif

   do j=Yed-1,Yed
   do k=Zst,Zed
   do i=0,Xmax+1
      buffer_sdy(i,k,j-Yed+2)=var(i,j,k)
   enddo
   enddo
   enddo
   call MPI_Sendrecv(buffer_sdy(0,Zst,1),(Xmax+2)*(Zed-Zst+1)*2,MPI_DOUBLE_PRECISION,src,8000,&
   &    buffer_rvy(0,Zst,1),(Xmax+2)*(Zed-Zst+1)*2,MPI_DOUBLE_PRECISION,dst,8000,MPI_COMM_WORLD,stat,ierr)
  if(south>=0)then
   do j=Yst-2,Yst-1
   do k=Zst,Zed
   do i=0,Xmax+1
      var(i,j,k)=buffer_rvy(i,k,j-Yst+3)
   enddo
   enddo
   enddo
  endif

   if(down>=0)then
        dst=down+idy*npz
   else
        dst=MPI_PROC_NULL
   endif
   if(top<npz)then
        src=top+idy*npz
   else
        src=MPI_PROC_NULL
   endif
   call MPI_Sendrecv(var(0,Yst-2,Zst),(Xmax+2)*(Yed-Yst+5)*2,MPI_DOUBLE_PRECISION,dst,9000,&
   &    var(0,Yst-2,Zed+1),(Xmax+2)*(Yed-Yst+5)*2,MPI_DOUBLE_PRECISION,src,9000,MPI_COMM_WORLD,stat,ierr)
   call MPI_Sendrecv(var(0,Yst-2,Zed-1),(Xmax+2)*(Yed-Yst+5)*2,MPI_DOUBLE_PRECISION,src,8000,&
   &    var(0,Yst-2,Zst-2),(Xmax+2)*(Yed-Yst+5)*2,MPI_DOUBLE_PRECISION,dst,8000,MPI_COMM_WORLD,stat,ierr)

end subroutine update

subroutine update_coefficient(i,j,k,u_0,v_0,flag1,flag2,flag3)
   use fangzhu
   implicit none

   integer,intent(in):: flag1(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
   integer,intent(in):: flag2(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
   integer,intent(in):: flag3(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)

   integer,intent(in):: i,j,k
   real(8)::u_0(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
   real(8)::v_0(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
   real(8)::ae,aw,an,as,at,ad,af!,aee,aww,ann,ass
   real(8)::diffe,diffw,diffn,diffs,difft,diffd
   real(8)::conve,convw,convn,convs,convt,convd
   integer :: east,west,north,south,top,down
   real(8) pb
   real(8) nx,ny,nz,uc,vc,wc,temp
   real*8 udefecorr1,udefecorr2,vdefecorr1,vdefecorr2,wdefecorr1,wdefecorr2

        east=i+1
        west=i-1
        south=j-1
        north=j+1
        top=k+1
        down=k-1

        a_u(i,j,k)=Vp/dt
        a_v(i,j,k)=Vp/dt
        a_w(i,j,k)=Vp/dt
        b_u(i,j,k)=0.0
        b_v(i,j,k)=0.0
        b_w(i,j,k)=0.0

        if(flag1(west,j,k)==INTERIOR)then
           diffw=nu*Sf1/dx
           if(mf(i,j,k,2)>0)then
                convw=mf(i,j,k,2)
                aw=mf(i,j,k,2)+diffw
                a_u(i,j,k)=a_u(i,j,k)+diffw
                a_v(i,j,k)=a_v(i,j,k)+diffw
                a_w(i,j,k)=a_w(i,j,k)+diffw
           else
                convw=0.0
                aw=diffw
                a_u(i,j,k)=a_u(i,j,k)-mf(i,j,k,2)+diffw
                a_v(i,j,k)=a_v(i,j,k)-mf(i,j,k,2)+diffw
                a_w(i,j,k)=a_w(i,j,k)-mf(i,j,k,2)+diffw
           endif
        else if(flag1(west,j,k)==V_INLET)then
           pb = p(i,j,k) - gradpx(i,j,k)*0.5*dx
           pf(i,j,k,2)=pb

           diffw=nu*Sf1/(0.5*dx)
           convw=0.0
           aw=0.0
           a_u(i,j,k)=a_u(i,j,k)+diffw
           a_v(i,j,k)=a_v(i,j,k)+diffw
           a_w(i,j,k)=a_w(i,j,k)+diffw
           b_u(i,j,k)=b_u(i,j,k)+diffw*uf(i,j,k,2)
           b_v(i,j,k)=b_v(i,j,k)+diffw*vf(i,j,k,2)
           b_w(i,j,k)=b_w(i,j,k)+diffw*wf(i,j,k,2)
           b_u(i,j,k)=b_u(i,j,k)+mf(i,j,k,2)*uf(i,j,k,2)
           b_v(i,j,k)=b_v(i,j,k)+mf(i,j,k,2)*vf(i,j,k,2)
           b_w(i,j,k)=b_w(i,j,k)+mf(i,j,k,2)*wf(i,j,k,2)
        else if(flag1(west,j,k)==WALL)then
           pb = p(i,j,k) - gradpx(i,j,k)*0.5*dx
           pf(i,j,k,2)=pb

           diffw=nu*Sf1/(0.5*dx)
           convw=0.0
           aw=0.0

           nx=-1.0
           ny=0.0
           nz=0.0

           a_u(i,j,k)=a_u(i,j,k)+diffw*(1.0-nx*nx)
           a_v(i,j,k)=a_v(i,j,k)+diffw*(1.0-ny*ny)
           a_w(i,j,k)=a_w(i,j,k)+diffw*(1.0-nz*nz)
        else
           print*,'unkown west bdy!,i,j,k:',i,j,k
        endif

        if(flag1(east,j,k)==INTERIOR)then
          diffe=nu*Sf1/dx
          if(mf(i+1,j,k,2)>0)then
                conve=-mf(i+1,j,k,2)
                ae=diffe
                a_u(i,j,k)=a_u(i,j,k)+mf(i+1,j,k,2)+diffe
                a_v(i,j,k)=a_v(i,j,k)+mf(i+1,j,k,2)+diffe
                a_w(i,j,k)=a_w(i,j,k)+mf(i+1,j,k,2)+diffe
          else
                conve=0.0
                ae=diffe-mf(i+1,j,k,2)
                a_u(i,j,k)=a_u(i,j,k)+diffe
                a_v(i,j,k)=a_v(i,j,k)+diffe
                a_w(i,j,k)=a_w(i,j,k)+diffe
          endif
        else if(flag1(east,j,k)==V_INLET)then
          diffe=nu*Sf1/(0.5*dx)
          conve=0.0
          ae=0.0
          a_u(i,j,k)=a_u(i,j,k)+diffe
          a_v(i,j,k)=a_v(i,j,k)+diffe
          a_w(i,j,k)=a_w(i,j,k)+diffe
          pb = p(i,j,k) + gradpx(i,j,k)*0.5*dx
          pf(i+1,j,k,2)=pb
          b_u(i,j,k)=b_u(i,j,k)+diffe*uf(i+1,j,k,2)
          b_v(i,j,k)=b_v(i,j,k)+diffe*vf(i+1,j,k,2)
          b_w(i,j,k)=b_w(i,j,k)+diffe*wf(i+1,j,k,2)
          b_u(i,j,k)=b_u(i,j,k)-mf(i+1,j,k,2)*uf(i+1,j,k,2)
          b_v(i,j,k)=b_v(i,j,k)-mf(i+1,j,k,2)*vf(i+1,j,k,2)
          b_w(i,j,k)=b_w(i,j,k)-mf(i+1,j,k,2)*wf(i+1,j,k,2)
        else if(flag1(east,j,k)==WALL)then
          pb = p(i,j,k) + gradpx(i,j,k)*0.5*dx
          pf(i+1,j,k,2)=pb

          nx=1.0
          ny=0.0
          nz=0.0
          diffe=nu*Sf1/(0.5*dx)
          conve=0.0
          ae=0.0

          a_u(i,j,k)=a_u(i,j,k)+diffe*(1.0 - nx * nx)
          a_v(i,j,k)=a_v(i,j,k)+diffe*(1.0 - ny * ny)
          a_w(i,j,k)=a_w(i,j,k)+diffe*(1.0 - nz * nz)

        else
           print*,'unkown east bdy!,i,j,k:',i,j,k
        endif


        if(flag1(i,south,k)==INTERIOR)then
          diffs=nu*Sf2/dy
          if(mf(i,j,k,1)>0)then
                convs=mf(i,j,k,1)
                as=mf(i,j,k,1)+diffs
                a_u(i,j,k)=a_u(i,j,k)+diffs
                a_v(i,j,k)=a_v(i,j,k)+diffs
                a_w(i,j,k)=a_w(i,j,k)+diffs
          else
                convs=0.0
                as=diffs
                a_u(i,j,k)=a_u(i,j,k)-mf(i,j,k,1)+diffs
                a_v(i,j,k)=a_v(i,j,k)-mf(i,j,k,1)+diffs
                a_w(i,j,k)=a_w(i,j,k)-mf(i,j,k,1)+diffs
          endif
        else if(flag1(i,south,k)==V_INLET)then
          diffs=nu*Sf2/(0.5*dy)
          convs=0.0
          as=0.0
          a_u(i,j,k)=a_u(i,j,k)+diffs
          a_v(i,j,k)=a_v(i,j,k)+diffs
          a_w(i,j,k)=a_w(i,j,k)+diffs
          pb = p(i,j,k) - gradpy(i,j,k)*0.5*dy
          pf(i,j,k,1)=pb
          b_u(i,j,k)=b_u(i,j,k)+diffs*uf(i,j,k,1)
          b_v(i,j,k)=b_v(i,j,k)+diffs*vf(i,j,k,1)
          b_w(i,j,k)=b_w(i,j,k)+diffs*wf(i,j,k,1)
          b_u(i,j,k)=b_u(i,j,k)+mf(i,j,k,1)*uf(i,j,k,1)
          b_v(i,j,k)=b_v(i,j,k)+mf(i,j,k,1)*vf(i,j,k,1)
          b_w(i,j,k)=b_w(i,j,k)+mf(i,j,k,1)*wf(i,j,k,1)
        else if(flag1(i,south,k)==WALL)then
          diffs=nu*Sf2/(0.5*dy)
          convs=0.0
          as=0.0

          nx=0.0
          ny=-1.0
          nz=0.0

          a_u(i,j,k)=a_u(i,j,k)+diffs*(1.0 - nx * nx)
          a_v(i,j,k)=a_v(i,j,k)+diffs*(1.0 - ny * ny)
          a_w(i,j,k)=a_w(i,j,k)+diffs*(1.0 - nz * nz)
          pb = p(i,j,k) - gradpy(i,j,k)*0.5*dy
          pf(i,j,k,1)=pb

        else
           print*,'unkown south bdy!,i,j,k:',i,j,k
        endif
        
        if(flag1(i,north,k)==INTERIOR)then
          diffn=nu*Sf2/dy
          
          if(mf(i,j+1,k,1)>0)then
                convn=-mf(i,j+1,k,1)
                an=diffn
                a_u(i,j,k)=a_u(i,j,k)+mf(i,j+1,k,1)+diffn
                a_v(i,j,k)=a_v(i,j,k)+mf(i,j+1,k,1)+diffn
                a_w(i,j,k)=a_w(i,j,k)+mf(i,j+1,k,1)+diffn
          else
                convn=0.0
                an=diffn-mf(i,j+1,k,1)
                a_u(i,j,k)=a_u(i,j,k)+diffn
                a_v(i,j,k)=a_v(i,j,k)+diffn
                a_w(i,j,k)=a_w(i,j,k)+diffn
          endif
        else if(flag1(i,north,k)==V_INLET)then
          diffn=nu*Sf2/(0.5*dy)
          convn=0.0
          an=0.0
          a_u(i,j,k)=a_u(i,j,k)+diffn
          a_v(i,j,k)=a_v(i,j,k)+diffn
          a_w(i,j,k)=a_w(i,j,k)+diffn
          pb = p(i,j,k) + gradpy(i,j,k)*0.5*dy
          pf(i,j+1,k,1)=pb
          b_u(i,j,k)=b_u(i,j,k)+diffn*uf(i,j+1,k,1)
          b_v(i,j,k)=b_v(i,j,k)+diffn*vf(i,j+1,k,1)
          b_w(i,j,k)=b_w(i,j,k)+diffn*wf(i,j+1,k,1)
          b_u(i,j,k)=b_u(i,j,k)-mf(i,j+1,k,1)*uf(i,j+1,k,1)
          b_v(i,j,k)=b_v(i,j,k)-mf(i,j+1,k,1)*vf(i,j+1,k,1)
          b_w(i,j,k)=b_w(i,j,k)-mf(i,j+1,k,1)*wf(i,j+1,k,1)
        else if(flag1(i,north,k)==WALL)then
          diffn=nu*Sf2/(0.5*dy)
          convn=0.0
          an=0.0

          nx=0.0
          ny=1.0
          nz=0.0

          a_u(i,j,k)=a_u(i,j,k)+diffn*(1.0 - nx * nx)
          a_v(i,j,k)=a_v(i,j,k)+diffn*(1.0 - ny * ny)
          a_w(i,j,k)=a_w(i,j,k)+diffn*(1.0 - nz * nz)
          pb = p(i,j,k) + gradpy(i,j,k)*0.5*dy
          pf(i,j+1,k,1)=pb
        else
           print*,'unkown north bdy!,i,j,k:',i,j,k
        endif

        if(flag1(i,j,down)==INTERIOR)then
          diffd=nu*Sf3/dz
          if(mf(i,j,k,3)>0)then
                convd=mf(i,j,k,3)
                ad=mf(i,j,k,3)+diffd
                a_u(i,j,k)=a_u(i,j,k)+diffd
                a_v(i,j,k)=a_v(i,j,k)+diffd
                a_w(i,j,k)=a_w(i,j,k)+diffd
          else
                convd=0.0
                ad=diffd
                a_u(i,j,k)=a_u(i,j,k)-mf(i,j,k,3)+diffd
                a_v(i,j,k)=a_v(i,j,k)-mf(i,j,k,3)+diffd
                a_w(i,j,k)=a_w(i,j,k)-mf(i,j,k,3)+diffd
          endif
        else if(flag1(i,j,down)==V_INLET)then
          diffd=nu*Sf3/(0.5*dz)
          convd=0.0
          ad=0.0
          a_u(i,j,k)=a_u(i,j,k)+diffd
          a_v(i,j,k)=a_v(i,j,k)+diffd
          a_w(i,j,k)=a_w(i,j,k)+diffd
          pb = p(i,j,k) - gradpz(i,j,k)*0.5*dz
          pf(i,j,k,3)=pb
          b_u(i,j,k)=b_u(i,j,k)+diffd*uf(i,j,k,3)
          b_v(i,j,k)=b_v(i,j,k)+diffd*vf(i,j,k,3)
          b_w(i,j,k)=b_w(i,j,k)+diffd*wf(i,j,k,3)
          b_u(i,j,k)=b_u(i,j,k)+mf(i,j,k,3)*uf(i,j,k,3)
          b_v(i,j,k)=b_v(i,j,k)+mf(i,j,k,3)*vf(i,j,k,3)
          b_w(i,j,k)=b_w(i,j,k)+mf(i,j,k,3)*wf(i,j,k,3)
        else if(flag1(i,j,down)==WALL)then
          diffd=nu*Sf3/(0.5*dz)
          convd=0.0
          ad=0.0

          nx=0.0
          ny=0.0
          nz=-1.0

          a_u(i,j,k)=a_u(i,j,k)+diffd*(1.0 - nx * nx)
          a_v(i,j,k)=a_v(i,j,k)+diffd*(1.0 - ny * ny)
          a_w(i,j,k)=a_w(i,j,k)+diffd*(1.0 - nz * nz)
          pb = p(i,j,k) - gradpz(i,j,k)*0.5*dz
          pf(i,j,k,3)=pb
        else
           print*,'unkown down bdy!,i,j,k:',i,j,k
        endif
        
        if(flag1(i,j,top)==INTERIOR)then
          difft=nu*Sf3/dz
          if(mf(i,j,k+1,3)>0)then
                convt=-mf(i,j,k+1,3)
                at=difft
                a_u(i,j,k)=a_u(i,j,k)+mf(i,j,k+1,3)+difft
                a_v(i,j,k)=a_v(i,j,k)+mf(i,j,k+1,3)+difft
                a_w(i,j,k)=a_w(i,j,k)+mf(i,j,k+1,3)+difft
          else
                convt=0.0
                at=difft-mf(i,j,k+1,3)
                a_u(i,j,k)=a_u(i,j,k)+difft
                a_v(i,j,k)=a_v(i,j,k)+difft
                a_w(i,j,k)=a_w(i,j,k)+difft
          endif
        else if(flag1(i,j,top)==V_INLET)then
          difft=nu*Sf3/(0.5*dz)
          convt=0.0
          at=0.0
          a_u(i,j,k)=a_u(i,j,k)+difft
          a_v(i,j,k)=a_v(i,j,k)+difft
          a_w(i,j,k)=a_w(i,j,k)+difft
          pb = p(i,j,k) + gradpz(i,j,k)*0.5*dz
          pf(i,j,k+1,3)=pb
          b_u(i,j,k)=b_u(i,j,k)+difft*uf(i,j,k+1,3)
          b_v(i,j,k)=b_v(i,j,k)+difft*vf(i,j,k+1,3)
          b_w(i,j,k)=b_w(i,j,k)+difft*wf(i,j,k+1,3)
          b_u(i,j,k)=b_u(i,j,k)-mf(i,j,k+1,3)*uf(i,j,k+1,3)
          b_v(i,j,k)=b_v(i,j,k)-mf(i,j,k+1,3)*vf(i,j,k+1,3)
          b_w(i,j,k)=b_w(i,j,k)-mf(i,j,k+1,3)*wf(i,j,k+1,3)
!        print*,'myid,i,j,k,top,b_v:',myid,i,j,k,top,b_v(i,j,k)
        else if(flag1(i,j,top)==WALL)then
          difft=nu*Sf3/(0.5*dz)
          convt=0.0
          at=0.0
           nx=0.0
           ny=0.0
           nz=1.0
          a_u(i,j,k)=a_u(i,j,k)+difft*(1.0 - nx * nx)
          a_v(i,j,k)=a_v(i,j,k)+difft*(1.0 - ny * ny)
          a_w(i,j,k)=a_w(i,j,k)+difft*(1.0 - nz * nz)
          pb = p(i,j,k) + gradpz(i,j,k)*0.5*dz
          pf(i,j,k+1,3)=pb
        else
           print*,'unkown north bdy!,i,j,k:',i,j,k
        endif

        b_u(i,j,k)=b_u(i,j,k)-Vp*(pf(i+1,j,k,2)-pf(i,j,k,2))/dx+Vp/dt*ui(i,j,k)
        b_v(i,j,k)=b_v(i,j,k)-Vp*(pf(i,j+1,k,1)-pf(i,j,k,1))/dy+Vp/dt*vi(i,j,k)
        b_w(i,j,k)=b_w(i,j,k)-Vp*(pf(i,j,k+1,3)-pf(i,j,k,3))/dz+Vp/dt*wi(i,j,k)

        a_uw(i,j,k)=aw
        a_ue(i,j,k)=ae
        a_un(i,j,k)=an
        a_us(i,j,k)=as
        a_ud(i,j,k)=ad
        a_ut(i,j,k)=at

        a_vw(i,j,k)=aw
        a_ve(i,j,k)=ae
        a_vn(i,j,k)=an
        a_vs(i,j,k)=as
        a_vd(i,j,k)=ad
        a_vt(i,j,k)=at

        a_ww(i,j,k)=aw
        a_we(i,j,k)=ae
        a_wn(i,j,k)=an
        a_ws(i,j,k)=as
        a_wd(i,j,k)=ad
        a_wt(i,j,k)=at

end subroutine update_coefficient
!bak 
 subroutine sub1_momentum(flag1,flag2,flag3)
    use fangzhu
    implicit none
 
    integer flag1(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
    integer flag2(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
    integer flag3(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
 
    integer :: i,j,k,l,it
    real*8 resi,err
    real(8)::u1,v1,w1!,D1,D2
    real(8) Df
    real(8)::gradpf1,gradpf2
    real(8)::a0,a1,pf_0,pf_1
    integer east,west,north,south,top,down
    integer numb
    integer::api(6)
    real(8)::apn_u(6)
    real(8)::apn_v(6)
    real(8)::apn_w(6)
    integer elem

!        if(myid==0) print*,'11111111111111'
        call flush(6)
        call MPI_Barrier(MPI_COMM_WORLD,ierr)

 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
         east=i+1
         west=i-1
        if(flag1(i,j,k)==INTERIOR)then
         if((flag1(east,j,k).eq.INTERIOR).and.(flag1(west,j,k).eq.INTERIOR))then
            gradpx(i,j,k)=(p(i+1,j,k)-p(i-1,j,k))/dx*0.5
            gradux(i,j,k)=(u(i+1,j,k)-u(i-1,j,k))/dx*0.5
            gradvx(i,j,k)=(v(i+1,j,k)-v(i-1,j,k))/dx*0.5
            gradwx(i,j,k)=(w(i+1,j,k)-w(i-1,j,k))/dx*0.5
         else if(flag1(east,j,k).ne.INTERIOR)then
            gradux(i,j,k)=(uf(i+1,j,k,2)-uf(i,j,k,2))/dx
            gradvx(i,j,k)=(vf(i+1,j,k,2)-vf(i,j,k,2))/dx
            gradwx(i,j,k)=(wf(i+1,j,k,2)-wf(i,j,k,2))/dx
            gradpx(i,j,k)=(pf(i+1,j,k,2)-pf(i,j,k,2))/dx
         else
            gradux(i,j,k)=(uf(i+1,j,k,2)-uf(i,j,k,2))/dx
            gradvx(i,j,k)=(vf(i+1,j,k,2)-vf(i,j,k,2))/dx
            gradwx(i,j,k)=(wf(i+1,j,k,2)-wf(i,j,k,2))/dx
            gradpx(i,j,k)=(pf(i+1,j,k,2)-pf(i,j,k,2))/dx
         endif
        endif
      enddo
      enddo
      enddo
 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
 
         north=j+1
         south=j-1
        if(flag1(i,j,k)==INTERIOR)then
         if((flag1(i,north,k).eq.INTERIOR).and.(flag1(i,south,k).eq.INTERIOR))then
            gradpy(i,j,k)=(p(i,j+1,k)-p(i,j-1,k))/dy*0.5
            graduy(i,j,k)=(u(i,j+1,k)-u(i,j-1,k))/dy*0.5
            gradvy(i,j,k)=(v(i,j+1,k)-v(i,j-1,k))/dy*0.5
            gradwy(i,j,k)=(w(i,j+1,k)-w(i,j-1,k))/dy*0.5
         else if(flag1(i,north,k).ne.INTERIOR)then
            graduy(i,j,k)=(uf(i,j+1,k,1)-uf(i,j,k,1))/dy
            gradvy(i,j,k)=(vf(i,j+1,k,1)-vf(i,j,k,1))/dy
            gradwy(i,j,k)=(wf(i,j+1,k,1)-wf(i,j,k,1))/dy
            gradpy(i,j,k)=(pf(i,j+1,k,1)-pf(i,j,k,1))/dy
         else
            graduy(i,j,k)=(uf(i,j+1,k,1)-uf(i,j,k,1))/dy
            gradvy(i,j,k)=(vf(i,j+1,k,1)-vf(i,j,k,1))/dy
            gradwy(i,j,k)=(wf(i,j+1,k,1)-wf(i,j,k,1))/dy
            gradpy(i,j,k)=(pf(i,j+1,k,1)-pf(i,j,k,1))/dy
         endif
        endif
      enddo
      enddo
      enddo
 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
 
         top=k+1
         down=k-1
        if(flag1(i,j,k)==INTERIOR)then
         if((flag1(i,j,top).eq.INTERIOR).and.(flag1(i,j,down).eq.INTERIOR))then
            gradpz(i,j,k)=(p(i,j,k+1)-p(i,j,k-1))/dz*0.5
            graduz(i,j,k)=(u(i,j,k+1)-u(i,j,k-1))/dz*0.5
            gradvz(i,j,k)=(v(i,j,k+1)-v(i,j,k-1))/dz*0.5
            gradwz(i,j,k)=(w(i,j,k+1)-w(i,j,k-1))/dz*0.5
         else if(flag1(i,j,top).ne.INTERIOR)then
            graduz(i,j,k)=(uf(i,j,k+1,3)-uf(i,j,k,3))/dz
            gradvz(i,j,k)=(vf(i,j,k+1,3)-vf(i,j,k,3))/dz
            gradwz(i,j,k)=(wf(i,j,k+1,3)-wf(i,j,k,3))/dz
            gradpz(i,j,k)=(pf(i,j,k+1,3)-pf(i,j,k,3))/dz
         else
            graduz(i,j,k)=(uf(i,j,k+1,3)-uf(i,j,k,3))/dz
            gradvz(i,j,k)=(vf(i,j,k+1,3)-vf(i,j,k,3))/dz
            gradwz(i,j,k)=(wf(i,j,k+1,3)-wf(i,j,k,3))/dz
            gradpz(i,j,k)=(pf(i,j,k+1,3)-pf(i,j,k,3))/dz
         endif
        endif
      enddo
      enddo
      enddo

      call update(gradux)
      call update(gradvx)
      call update(gradwx)
      call update(gradpx)


      call update(graduy)
      call update(gradvy)
      call update(gradwy)
      call update(gradpy)


      call update(graduz)
      call update(gradvz)
      call update(gradwz)
      call update(gradpz)

    b_u=0.0
    a_u=0.0
    a_uw=0.0
    a_ue=0.0
    a_us=0.0
    a_un=0.0
    a_ut=0.0
    a_ud=0.0
 
    b_v=0.0
    a_v=0.0
    a_vw=0.0
    a_ve=0.0
    a_vs=0.0
    a_vn=0.0
    a_vt=0.0
    a_vd=0.0
 
    b_w=0.0
    a_w=0.0
    a_ww=0.0
    a_we=0.0
    a_ws=0.0
    a_wn=0.0
    a_wt=0.0
    a_wd=0.0
 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
        if(flag1(i,j,k)==INTERIOR)then
         call update_coefficient(i,j,k,u,v,flag1,flag2,flag3)
        endif
      enddo
      enddo
      enddo

      call update(a_u)
      call update(a_v)
      call update(a_w)

      call mem_alloca_u(gridsnum)
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
        numb=0
        apn_u=0.0
        apn_v=0.0
        apn_w=0.0
        api=0
        if(flag1(i,j,k)==INTERIOR)then
          if(flag1(i-1,j,k)==INTERIOR)then
             numb=numb+1
             api(numb)=(i-1)+(j-1)*Xmax+(k-1)*Ymax*Xmax
             apn_u(numb)=-a_uw(i,j,k)
             apn_v(numb)=-a_vw(i,j,k)
             apn_w(numb)=-a_ww(i,j,k)
          endif
          if(flag1(i+1,j,k)==INTERIOR)then
             numb=numb+1
             api(numb)=(i+1)+(j-1)*Xmax+(k-1)*Ymax*Xmax
             apn_u(numb)=-a_ue(i,j,k)
             apn_v(numb)=-a_ve(i,j,k)
             apn_w(numb)=-a_we(i,j,k)
          endif
          if(flag1(i,j-1,k)==INTERIOR)then
             numb=numb+1
             api(numb)=i+(j-2)*Xmax+(k-1)*Ymax*Xmax
             apn_u(numb)=-a_us(i,j,k)
             apn_v(numb)=-a_vs(i,j,k)
             apn_w(numb)=-a_ws(i,j,k)
          endif
          if(flag1(i,j+1,k)==INTERIOR)then
             numb=numb+1
             api(numb)=i+j*Xmax+(k-1)*Ymax*Xmax
             apn_u(numb)=-a_un(i,j,k)
             apn_v(numb)=-a_vn(i,j,k)
             apn_w(numb)=-a_wn(i,j,k)
          endif
          if(flag1(i,j,k-1)==INTERIOR)then
             numb=numb+1
             api(numb)=i+(j-1)*Xmax+(k-2)*Ymax*Xmax
             apn_u(numb)=-a_ud(i,j,k)
             apn_v(numb)=-a_vd(i,j,k)
             apn_w(numb)=-a_wd(i,j,k)
          endif
          if(flag1(i,j,k+1)==INTERIOR)then
             numb=numb+1
             api(numb)=i+(j-1)*Xmax+k*Ymax*Xmax
             apn_u(numb)=-a_ut(i,j,k)
             apn_v(numb)=-a_vt(i,j,k)
             apn_w(numb)=-a_wt(i,j,k)
          endif
        endif
        elem=i+(j-1)*Xmax+(k-1)*Ymax*Xmax
        call Q_Set_U(elem,numb,a_u(i,j,k),a_v(i,j,k),a_w(i,j,k),b_u(i,j,k),b_v(i,j,k),b_w(i,j,k),api,apn_u,apn_v,apn_w)
      enddo
      enddo
      enddo


     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
        elem=i+(j-1)*Xmax+(k-1)*Ymax*Xmax
        call SETUVW(u(i,j,k),v(i,j,k),w(i,j,k),elem)
     enddo
     enddo
     enddo
        
      call SOLVEMatrix_fasp_U()
     if (mod(istep,NUMAX)==0) then
        call output_uvw(istep)
     endif
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
        elem=i+(j-1)*Xmax+(k-1)*Ymax*Xmax
        call GETUVW(u(i,j,k),v(i,j,k),w(i,j,k),elem)
     enddo
     enddo
     enddo

      call mem_dealloca_u()
 
      call update(u)
 
      call update(v)

      call update(w)
 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
         east=i+1
         west=i-1
        if(flag1(i,j,k)==INTERIOR)then
         if((flag1(east,j,k).eq.INTERIOR).and.(flag1(west,j,k).eq.INTERIOR))then
            gradpx(i,j,k)=(p(i+1,j,k)-p(i-1,j,k))/dx*0.5
         else if(flag1(east,j,k).ne.INTERIOR)then
            gradpx(i,j,k)=(pf(i+1,j,k,2)-pf(i,j,k,2))/dx
         else
            gradpx(i,j,k)=(pf(i+1,j,k,2)-pf(i,j,k,2))/dx
         endif
        endif
      enddo
      enddo
      enddo
      call update(gradpx)
 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
         north=j+1
         south=j-1
        if(flag1(i,j,k)==INTERIOR)then
         if((flag1(i,north,k).eq.INTERIOR).and.(flag1(i,south,k).eq.INTERIOR))then
            gradpy(i,j,k)=(p(i,j+1,k)-p(i,j-1,k))/dy*0.5
         else if(flag1(i,north,k).ne.INTERIOR)then
            gradpy(i,j,k)=(pf(i,j+1,k,1)-pf(i,j,k,1))/dy
         else
            gradpy(i,j,k)=(pf(i,j+1,k,1)-pf(i,j,k,1))/dy
         endif
        endif
      enddo
      enddo
      enddo
      call update(gradpy)
 
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
 
         top=k+1
         down=k-1
        if(flag1(i,j,k)==INTERIOR)then
         if((flag1(i,j,top).eq.INTERIOR).and.(flag1(i,j,down).eq.INTERIOR))then
            gradpz(i,j,k)=(p(i,j,k+1)-p(i,j,k-1))/dz*0.5
         else if(flag1(i,j,top).ne.INTERIOR)then
            gradpz(i,j,k)=(pf(i,j,k+1,3)-pf(i,j,k,3))/dz
         else
            gradpz(i,j,k)=(pf(i,j,k+1,3)-pf(i,j,k,3))/dz
         endif
        endif
      enddo
      enddo
      enddo
      call update(gradpz)
 
 
  l=1
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if(flag1(i,j,k)==INTERIOR)then
         if(flag1(i,j-1,k)==INTERIOR)then
            u1=0.5*(u(i,j-1,k)+u(i,j,k))
            v1=0.5*(v(i,j-1,k)+v(i,j,k))
            w1=0.5*(w(i,j-1,k)+w(i,j,k))
            Df=2*Vp/(a_v(i,j-1,k)+a_v(i,j,k))
            gradpf1=(p(i,j,k)-p(i,j-1,k))/dy
            gradpf2=0.5*(gradpy(i,j-1,k)+gradpy(i,j,k))
            vf(i,j,k,l)=v1-Df*(gradpf1-gradpf2)
            mf(i,j,k,l)=vf(i,j,k,l)*Sf2
            uf(i,j,k,l)=u1
            wf(i,j,k,l)=w1
         endif
         if(flag1(i,j+1,k)==INTERIOR)then
            u1=0.5*(u(i,j,k)+u(i,j+1,k))
            v1=0.5*(v(i,j,k)+v(i,j+1,k))
            w1=0.5*(w(i,j,k)+w(i,j+1,k))
            Df=2*Vp/(a_v(i,j,k)+a_v(i,j+1,k))
            gradpf1=(p(i,j+1,k)-p(i,j,k))/dy
            gradpf2=0.5*(gradpy(i,j,k)+gradpy(i,j+1,k))
            vf(i,j+1,k,l)=v1-Df*(gradpf1-gradpf2)
            mf(i,j+1,k,l)=vf(i,j+1,k,l)*Sf2
            uf(i,j+1,k,l)=u1
            wf(i,j+1,k,l)=w1
         endif
       endif
     enddo
     enddo
     enddo
 
  l=2
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if(flag1(i,j,k)==INTERIOR)then
         if(flag1(i-1,j,k)==INTERIOR)then
            u1=0.5*(u(i-1,j,k)+u(i,j,k))
            v1=0.5*(v(i-1,j,k)+v(i,j,k))
            w1=0.5*(w(i-1,j,k)+w(i,j,k))
            Df=2*Vp/(a_u(i-1,j,k)+a_u(i,j,k))
            gradpf1=(p(i,j,k)-p(i-1,j,k))/dx
            gradpf2=0.5*(gradpx(i-1,j,k)+gradpx(i,j,k))
            uf(i,j,k,l)=u1-Df*(gradpf1-gradpf2)
            mf(i,j,k,l)=uf(i,j,k,l)*Sf1
            vf(i,j,k,l)=v1
            wf(i,j,k,l)=w1
         endif
         if(flag1(i+1,j,k)==INTERIOR)then
            u1=0.5*(u(i,j,k)+u(i+1,j,k))
            v1=0.5*(v(i,j,k)+v(i+1,j,k))
            w1=0.5*(w(i,j,k)+w(i+1,j,k))
            Df=2*Vp/(a_u(i,j,k)+a_u(i+1,j,k))
            gradpf1=(p(i+1,j,k)-p(i,j,k))/dx
            gradpf2=0.5*(gradpx(i,j,k)+gradpx(i+1,j,k))
            uf(i+1,j,k,l)=u1-Df*(gradpf1-gradpf2)
            mf(i+1,j,k,l)=uf(i+1,j,k,l)*Sf1
            vf(i+1,j,k,l)=v1
            wf(i+1,j,k,l)=w1
         endif
       endif
     enddo
     enddo
     enddo

  l=3
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if(flag1(i,j,k)==INTERIOR)then
         if(flag1(i,j,k-1)==INTERIOR)then
            u1=0.5*(u(i,j,k-1)+u(i,j,k))
            v1=0.5*(v(i,j,k-1)+v(i,j,k))
            w1=0.5*(w(i,j,k-1)+w(i,j,k))
            Df=2*Vp/(a_u(i,j,k-1)+a_u(i,j,k))
            gradpf1=(p(i,j,k)-p(i,j,k-1))/dz
            gradpf2=0.5*(gradpz(i,j,k-1)+gradpz(i,j,k))
            wf(i,j,k,l)=w1-Df*(gradpf1-gradpf2)
            mf(i,j,k,l)=wf(i,j,k,l)*Sf3
            uf(i,j,k,l)=u1
            vf(i,j,k,l)=v1
         endif
         if(flag1(i,j,k+1)==INTERIOR)then
            u1=0.5*(u(i,j,k)+u(i,j,k+1))
            v1=0.5*(v(i,j,k)+v(i,j,k+1))
            w1=0.5*(w(i,j,k)+w(i,j,k+1))
            Df=2*Vp/(a_u(i,j,k)+a_u(i,j,k+1))
            gradpf1=(p(i,j,k+1)-p(i,j,k))/dz
            gradpf2=0.5*(gradpz(i,j,k)+gradpz(i,j,k+1))
            wf(i,j,k+1,l)=w1-Df*(gradpf1-gradpf2)
            mf(i,j,k+1,l)=wf(i,j,k+1,l)*Sf3
            uf(i,j,k+1,l)=u1
            vf(i,j,k+1,l)=v1
         endif
       endif
     enddo
     enddo
     enddo
 
 end subroutine sub1_momentum
 
 subroutine sub3_pressure(flag1,flag2,flag3)
    use fangzhu
    implicit none
 
    integer flag1(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
    integer flag2(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
    integer flag3(0:Xmax+1,Yst-2:Yed+2,Zst-2:Zed+2)
 
    integer i,j,k,l
    real(8) mfe,mfw,mfn,mfs,mfd,mft
    real(8) ufe,ufw,vfs,vfn,wfd,wft
    real(8) Dfe,Dfw,Dfn,Dfs,Dfd,Dft
    real(8) gradpe,gradpw,gradps,gradpn,gradpc,gradpd,gradpt
    real(8) gradpfe,gradpfw,gradpfs,gradpfn,gradpfd,gradpft
    real(8) gradpfe2,gradpfw2,gradpfs2,gradpfn2,gradpfd2,gradpft2
 
    real(8) error,sump0,sump,err
    real(8)::D1,D2
    real(8)::gradpcf1,gradpcf2
    real(8)::a0,a1,pf_0,pf_1
    integer east,west,north,south,down,top
    real(8) nx,ny,nz,uc,vc,wc
    integer numb
    integer::api(6)
    real(8)::apn_p(6)
    integer elem

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      call mem_alloca_p(gridsnum)

      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
        numb=0
        apn_p=0.0
        api=0

        if(flag2(i,j,k)==INTERIOR)then
         if(flag2(i-1,j,k)==INTERIOR)then
           Dfw=2*Vp/(a_u(i-1,j,k)+a_u(i,j,k))*Sf1/dx!0.5*(Dw+Dc)*Sf/dx
           mfw=mf(i,j,k,2)*(-1.0)

           numb=numb+1
           api(numb)=(i-1)+(j-1)*Xmax+(k-1)*Ymax*Xmax
           apn_p(numb)=-Dfw
         else
           Dfw=0.0
           mfw=mf(i,j,k,2)*(-1.0)
         endif
 
         if(flag2(i+1,j,k)==INTERIOR)then
           Dfe=2*Vp/(a_u(i+1,j,k)+a_u(i,j,k))*Sf1/dx!0.5*(De+Dc)*Sf/dx
           mfe=mf(i+1,j,k,2)

           numb=numb+1
           api(numb)=(i+1)+(j-1)*Xmax+(k-1)*Ymax*Xmax
           apn_p(numb)=-Dfe

         else
            mfe=mf(i+1,j,k,2)!0.0
            Dfe=0.0
         endif
 
         if(flag2(i,j-1,k)==INTERIOR)then
           Dfs=2*Vp/(a_v(i,j-1,k)+a_v(i,j,k))*Sf2/dy!0.5*(Ds+Dc)*Sf/dy
           mfs=mf(i,j,k,1)*(-1.0)

           numb=numb+1
           api(numb)=(i)+(j-2)*Xmax+(k-1)*Ymax*Xmax
           apn_p(numb)=-Dfs

         else
           Dfs=0.0
           mfs=mf(i,j,k,1)*(-1.0)!0.0
         endif
 
         if(flag2(i,j+1,k)==INTERIOR)then
           Dfn=2*Vp/(a_v(i,j+1,k)+a_v(i,j,k))*Sf2/dy!0.5*(Dn+Dc)*Sf/dy
           mfn=mf(i,j+1,k,1)*1.0

           numb=numb+1
           api(numb)=(i)+j*Xmax+(k-1)*Ymax*Xmax
           apn_p(numb)=-Dfn

         else
           Dfn=0.0
           mfn=mf(i,j+1,k,1)!0.0
         endif
 
         if(flag2(i,j,k-1)==INTERIOR)then
           Dfd=2*Vp/(a_w(i,j,k-1)+a_w(i,j,k))*Sf3/dz!0.5*(Ds+Dc)*Sf/dy
           mfd=mf(i,j,k,3)*(-1.0)

           numb=numb+1
           api(numb)=(i)+(j-1)*Xmax+(k-2)*Ymax*Xmax
           apn_p(numb)=-Dfd

         else
           Dfd=0.0
           mfd=mf(i,j,k,3)*(-1.0)!0.0
         endif
 
         if(flag2(i,j,k+1)==INTERIOR)then
           Dft=2*Vp/(a_w(i,j,k+1)+a_w(i,j,k))*Sf3/dz!0.5*(Dn+Dc)*Sf/dy
           mft=mf(i,j,k+1,3)*1.0

           numb=numb+1
           api(numb)=(i)+(j-1)*Xmax+k*Ymax*Xmax
           apn_p(numb)=-Dft

         else
           Dft=0.0
           mft=mf(i,j,k+1,3)!0.0
         endif

         b_p(i,j,k)=-(mfw+mfe+mfs+mfn+mfd+mft)
         a_p(i,j,k)=a_p(i,j,k)+(Dfw+Dfe+Dfs+Dfn+Dfd+Dft)
         if(abs(a_p(i,j,k))<1.0e-15)then
           print*,'Why???,i,j,k,a_p,Dfw,Dfe,Dfs,Dfn:',i,j,k,a_p(i,j,k),Dfw,Dfe,Dfs,Dfn,Dfd,Dft
         endif
       endif
        
        elem=i+(j-1)*Xmax+(k-1)*Ymax*Xmax
        call Q_Set_P(elem,numb,a_p(i,j,k),b_p(i,j,k),api,apn_p)

      enddo !i
      enddo !j
      enddo !k
 
      call SOLVEMatrix_fasp_P()
     if (mod(istep,NUMAX)==0) then
        call output_p(istep)
     endif
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
        elem=i+(j-1)*Xmax+(k-1)*Ymax*Xmax
        call getpc(pc(i,j,k),elem)
     enddo
     enddo
     enddo
     call mem_dealloca_p()
      call update(pc)

 
 !CorrectPressureFace
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
        if(flag2(i,j,k)==INTERIOR)then
 
         if(flag2(i-1,j,k)==V_INLET)then
           pcf(i,j,k,2)=pc(i,j,k)
         else if(flag2(i-1,j,k)==WALL)then
           pcf(i,j,k,2)=pc(i,j,k)
         else
         endif !flag1(i-1,j)
 
         if(flag2(i+1,j,k)==V_INLET)then
           pcf(i+1,j,k,2)=pc(i,j,k)
         else if(flag2(i+1,j,k)==WALL)then
           pcf(i+1,j,k,2)=pc(i,j,k)
         else
         endif !flag1(i+1,j)
 
         if(flag2(i,j-1,k)==V_INLET)then
           pcf(i,j,k,1)=pc(i,j,k)
         else if(flag2(i,j-1,k)==WALL)then
           pcf(i,j,k,1)=pc(i,j,k)
         else
         endif !flag1(i,j-1)
 
         if(flag2(i,j+1,k)==V_INLET)then
           pcf(i,j+1,k,1)=pc(i,j,k)
         else if(flag2(i,j+1,k)==WALL)then
           pcf(i,j+1,k,1)=pc(i,j,k)
         else
         endif !flag1(i,j+1)

         if(flag2(i,j,k-1)==V_INLET)then
           pcf(i,j,k,3)=pc(i,j,k)
         else if(flag2(i,j,k-1)==WALL)then
           pcf(i,j,k,3)=pc(i,j,k)
         else
         endif !flag1(i,j-1)
 
         if(flag2(i,j,k+1)==V_INLET)then
           pcf(i,j,k+1,3)=pc(i,j,k)
         else if(flag2(i,j,k+1)==WALL)then
           pcf(i,j,k+1,3)=pc(i,j,k)
         else
         endif !flag1(i,j+1)

        endif !flag1(i,j)==INTERIOR
     enddo
     enddo
     enddo
 
 !CorrectPressureField
     !对速度p做修正下一时刻的un大小
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
        if(flag3(i,j,k)==INTERIOR)then
          p(i,j,k)=p(i,j,k)+relaxp*pc(i,j,k)
        endif
     enddo
     enddo
     enddo
      call update(p)
!        call timer(ed)
!      if(myid==0) print*,'time new p:',ed-st
!        call timer(st)
 
 !面上的压强不是用0.5计算的
  l=1
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if((flag2(i,j,k)==INTERIOR.and.flag2(i,j-1,k)==INTERIOR))then
         a0=(a_u(i,j,k)+a_v(i,j,k)+a_w(i,j,k))/3.0
         a1=(a_u(i,j-1,k)+a_v(i,j-1,k)+a_w(i,j-1,k))/3.0
         pf_0=p(i,j,k)-gradpy(i,j,k)*0.5*dy
         pf_1=p(i,j-1,k)+gradpy(i,j-1,k)*0.5*dy
         pf(i,j,k,l)=(a0*pf_0+a1*pf_1)/(a0+a1)
       endif
       if((flag2(i,j+1,k)==INTERIOR.and.flag2(i,j,k)==INTERIOR))then
         a0=(a_u(i,j+1,k)+a_v(i,j+1,k)+a_w(i,j+1,k))/3.0
         a1=(a_u(i,j,k)+a_v(i,j,k)+a_w(i,j,k))/3.0
         pf_0=p(i,j+1,k)-gradpy(i,j+1,k)*0.5*dy
         pf_1=p(i,j,k)+gradpy(i,j,k)*0.5*dy
         pf(i,j+1,k,l)=(a0*pf_0+a1*pf_1)/(a0+a1)
       endif
     enddo
     enddo
     enddo

  l=2
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if((flag2(i,j,k)==INTERIOR.and.flag2(i-1,j,k)==INTERIOR))then
         a0=(a_u(i,j,k)+a_v(i,j,k)+a_w(i,j,k))/3.0
         a1=(a_u(i-1,j,k)+a_v(i-1,j,k)+a_w(i-1,j,k))/3.0
         pf_0=p(i,j,k)-gradpx(i,j,k)*0.5*dx
         pf_1=p(i-1,j,k)+gradpx(i-1,j,k)*0.5*dx
         pf(i,j,k,l)=(a0*pf_0+a1*pf_1)/(a0+a1)
       endif
       if((flag2(i+1,j,k)==INTERIOR.and.flag2(i,j,k)==INTERIOR))then
         a0=(a_u(i+1,j,k)+a_v(i+1,j,k)+a_w(i+1,j,k))/3.0
         a1=(a_u(i,j,k)+a_v(i,j,k)+a_w(i,j,k))/3.0
         pf_0=p(i+1,j,k)-gradpx(i+1,j,k)*0.5*dx
         pf_1=p(i,j,k)+gradpx(i,j,k)*0.5*dx
         pf(i+1,j,k,l)=(a0*pf_0+a1*pf_1)/(a0+a1)
       endif
     enddo
     enddo
     enddo

  l=3
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if((flag2(i,j,k)==INTERIOR.and.flag2(i,j,k-1)==INTERIOR))then
         a0=(a_u(i,j,k)+a_v(i,j,k)+a_w(i,j,k))/3.0
         a1=(a_u(i,j,k-1)+a_v(i,j,k-1)+a_w(i,j,k-1))/3.0
         pf_0=p(i,j,k)-gradpz(i,j,k)*0.5*dz
         pf_1=p(i,j,k-1)+gradpz(i,j,k-1)*0.5*dz
         pf(i,j,k,l)=(a0*pf_0+a1*pf_1)/(a0+a1)
       endif
       if((flag2(i,j,k+1)==INTERIOR.and.flag2(i,j,k)==INTERIOR))then
         a0=(a_u(i,j,k+1)+a_v(i,j,k+1)+a_w(i,j,k+1))/3.0
         a1=(a_u(i,j,k)+a_v(i,j,k)+a_w(i,j,k))/3.0
         pf_0=p(i,j,k+1)-gradpz(i,j,k+1)*0.5*dz
         pf_1=p(i,j,k)+gradpz(i,j,k)*0.5*dz
         pf(i,j,k+1,l)=(a0*pf_0+a1*pf_1)/(a0+a1)
       endif
     enddo
     enddo
     enddo
!        call timer(ed)
!      if(myid==0) print*,'time new pf:',ed-st
!        call timer(st)
 
 !CorrectFaceFlux
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
         east=i+1
         west=i-1
         north=j+1
         south=j-1
         top=k+1
         down=k-1

        if(flag3(i,j,k)==INTERIOR)then
         if((flag2(east,j,k).eq.INTERIOR).and.(flag2(west,j,k).eq.INTERIOR))then
            gradpcx(i,j,k)=(pc(i+1,j,k)-pc(i-1,j,k))/dx*0.5
         else if(flag2(east,j,k).ne.INTERIOR)then
            gradpcx(i,j,k)=(pcf(i+1,j,k,2)-pcf(i,j,k,2))/dx
         else
            gradpcx(i,j,k)=(pcf(i+1,j,k,2)-pcf(i,j,k,2))/dx
         endif
 
         if((flag2(i,north,k).eq.INTERIOR).and.(flag2(i,south,k).eq.INTERIOR))then
            gradpcy(i,j,k)=(pc(i,j+1,k)-pc(i,j-1,k))/dy*0.5
         else if(flag2(i,north,k).ne.INTERIOR)then
            gradpcy(i,j,k)=(pcf(i,j+1,k,1)-pcf(i,j,k,1))/dy
         else
            gradpcy(i,j,k)=(pcf(i,j+1,k,1)-pcf(i,j,k,1))/dy
         endif
 
         if((flag2(i,j,top).eq.INTERIOR).and.(flag2(i,j,down).eq.INTERIOR))then
            gradpcz(i,j,k)=(pc(i,j,k+1)-pc(i,j,k-1))/dz*0.5
         else if(flag2(i,j,top).ne.INTERIOR)then
            gradpcz(i,j,k)=(pcf(i,j,k+1,3)-pcf(i,j,k,3))/dz
         else
            gradpcz(i,j,k)=(pcf(i,j,k+1,3)-pcf(i,j,k,3))/dz
         endif
        endif
      enddo
      enddo
      enddo

      call update(gradpcx)
      call update(gradpcy)
      call update(gradpcz)
!        call timer(ed)
!      if(myid==0) print*,'time gradpc:',ed-st
!        call timer(st)
 
      l=1
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
       if(flag3(i,j,k)==INTERIOR)then
         if(flag3(i,j-1,k)==INTERIOR)then
            gradpcf1=(pc(i,j,k)-pc(i,j-1,k))/dy
            vf(i,j,k,l)=vf(i,j,k,l)-2*Vp/(a_v(i,j-1,k)+a_v(i,j,k))*gradpcf1
            mf(i,j,k,l)=vf(i,j,k,l)*Sf2
         endif
       endif
      enddo
      enddo
      enddo
      j=Yed
      do k=Zst,Zed
      do i=1,Xmax
        if(flag3(i,j,k)==INTERIOR .and.flag3(i,j+1,k)==INTERIOR)then
!bak            gradpcf1=-(pc(i,j+1,k)-pc(i,j,k))/dy
            gradpcf1=(pc(i,j+1,k)-pc(i,j,k))/dy
            vf(i,j+1,k,l)=vf(i,j+1,k,l)-2*Vp/(a_v(i,j+1,k)+a_v(i,j,k))*gradpcf1
            mf(i,j+1,k,l)=vf(i,j+1,k,l)*Sf2
        endif
      enddo
      enddo
      
      l=2
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
       if(flag3(i,j,k)==INTERIOR)then
         if(flag3(i-1,j,k)==INTERIOR)then
            gradpcf1=(pc(i,j,k)-pc(i-1,j,k))/dx
            uf(i,j,k,l)=uf(i,j,k,l)-2.0*Vp/(a_u(i-1,j,k)+a_u(i,j,k))*gradpcf1
            mf(i,j,k,l)=uf(i,j,k,l)*Sf1
         endif
       endif
      enddo
      enddo
      enddo

      i=Xmax
      do k=Zst,Zed
      do j=Yst,Yed
        if((flag3(i,j,k)==INTERIOR).and.(flag3(i+1,j,k)==INTERIOR))then
            gradpcf1=(pc(i+1,j,k)-pc(i,j,k))/dx
            uf(i+1,j,k,l)=uf(i+1,j,k,l)-2.0*Vp/(a_u(i+1,j,k)+a_u(i,j,k))*gradpcf1
            mf(i+1,j,k,l)=uf(i+1,j,k,l)*Sf1
        endif
      enddo
      enddo
 
      l=3
      do k=Zst,Zed
      do j=Yst,Yed
      do i=1,Xmax
       if(flag3(i,j,k)==INTERIOR)then
         if(flag3(i,j,k-1)==INTERIOR)then
            gradpcf1=(pc(i,j,k)-pc(i,j,k-1))/dz
            wf(i,j,k,l)=wf(i,j,k,l)-2.0*Vp/(a_u(i,j,k-1)+a_u(i,j,k))*gradpcf1
            mf(i,j,k,l)=wf(i,j,k,l)*Sf3
         endif
       endif
      enddo
      enddo
      enddo

     k=Zed
     do j=Yst,Yed
     do i=1,Xmax
        if((flag3(i,j,k)==INTERIOR).and.(flag3(i,j,k+1)==INTERIOR))then
            gradpcf1=(pc(i,j,k+1)-pc(i,j,k))/dz
            wf(i,j,k+1,l)=wf(i,j,k+1,l)-2.0*Vp/(a_u(i,j,k+1)+a_u(i,j,k))*gradpcf1
            mf(i,j,k+1,l)=wf(i,j,k+1,l)*Sf3
        endif
     enddo
     enddo

!        call timer(ed)
!      if(myid==0) print*,'time new mf:',ed-st
!        call timer(st)
     !对速度u做修正下一时刻的un大小
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if(flag3(i,j,k)==INTERIOR)then
         u(i,j,k)=u(i,j,k)-gradpcx(i,j,k)*Vp/a_u(i,j,k)
       endif
     enddo
     enddo
     enddo

     call update(u)
 
     !对速度vn做修正下一时刻的un大小
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if(flag3(i,j,k)==INTERIOR)then
         v(i,j,k)=v(i,j,k)-gradpcy(i,j,k)*Vp/a_v(i,j,k)
       endif
     enddo
     enddo
     enddo
     call update(v)
 
     !对速度vn做修正下一时刻的un大小
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
       if(flag3(i,j,k)==INTERIOR)then
         w(i,j,k)=w(i,j,k)-gradpcz(i,j,k)*Vp/a_w(i,j,k)
       endif
     enddo
     enddo
     enddo
     call update(w)

 !CorrectFaceUVWP
  l=1
     do k=Zst,Zed
     do j=Yst,Yed
     do i=1,Xmax
        if(flag3(i,j,k)==INTERIOR)then
         if(flag3(i-1,j,k)==V_INLET)then
            pf(i,j,k,2)=pf(i,j,k,2)+pcf(i,j,k,2)*relaxp
         else if(flag3(i-1,j,k)==WALL)then
            mf(i,j,k,2)=0.0
            uf(i,j,k,2)=0.0
            vf(i,j,k,2)=0.0
            wf(i,j,k,2)=0.0
            pf(i,j,k,2)=pf(i,j,k,2)+pcf(i,j,k,2)*relaxp
         else
         endif !flag1(i-1,j)
 
         if(flag3(i+1,j,k)==V_INLET)then
            pf(i+1,j,k,2)=pf(i+1,j,k,2)+pcf(i+1,j,k,2)*relaxp
         else if(flag3(i+1,j,k)==WALL)then
            mf(i+1,j,k,2)=0.0
            uf(i+1,j,k,2)=0.0
            vf(i+1,j,k,2)=0.0
            wf(i+1,j,k,2)=0.0
            pf(i+1,j,k,2)=pf(i+1,j,k,2)+pcf(i+1,j,k,2)*relaxp
         else
         endif !flag1(i+1,j) +1
 
         if(flag3(i,j-1,k)==V_INLET)then
            pf(i,j,k,1)=pf(i,j,k,1)+pcf(i,j,k,1)*relaxp
         else if(flag3(i,j-1,k)==WALL)then
            mf(i,j,k,1)=0.0
            uf(i,j,k,1)=0.0
            vf(i,j,k,1)=0.0
            wf(i,j,k,1)=0.0
            pf(i,j,k,1)=pf(i,j,k,1)+pcf(i,j,k,1)*relaxp
         else
         endif !flag1(i,j-1)
 
         if(flag3(i,j+1,k)==V_INLET)then
            pf(i,j+1,k,1)=pf(i,j+1,k,1)+pcf(i,j+1,k,1)*relaxp
         else if(flag3(i,j+1,k)==WALL)then
            mf(i,j+1,k,1)=0.0
            uf(i,j+1,k,1)=0.0
            vf(i,j+1,k,1)=0.0
            wf(i,j+1,k,1)=0.0
            pf(i,j+1,k,1)=pf(i,j+1,k,1)+pcf(i,j+1,k,1)*relaxp
         else
         endif !flag1(i,j+1)
 
         if(flag3(i,j,k-1)==V_INLET)then
            pf(i,j,k,3)=pf(i,j,k,3)+pcf(i,j,k,3)*relaxp
         else if(flag3(i,j,k-1)==WALL)then
            mf(i,j,k,3)=0.0
            uf(i,j,k,3)=0.0
            vf(i,j,k,3)=0.0
            wf(i,j,k,3)=0.0
            pf(i,j,k,3)=pf(i,j,k,3)+pcf(i,j,k,3)*relaxp

         else
         endif !flag1(i,j-1)
 
         if(flag3(i,j,k+1)==V_INLET)then
            pf(i,j,k+1,3)=pf(i,j,k+1,3)+pcf(i,j,k+1,3)*relaxp
         else if(flag3(i,j,k+1)==WALL)then
            mf(i,j,k+1,3)=0.0
            uf(i,j,k+1,3)=0.0
            vf(i,j,k+1,3)=0.0
            wf(i,j,k+1,3)=0.0
            pf(i,j,k+1,3)=pf(i,j,k+1,3)+pcf(i,j,k+1,3)*relaxp
         else
         endif !flag1(i,j+1)

        endif !flag1(i,j)==INTERIOR
     enddo
     enddo
     enddo
!        call timer(ed)
!      if(myid==0) print*,'time new uv:',ed-st
 
 end subroutine sub3_pressure

subroutine dumpp()
   use fangzhu
   implicit none

   character*6 rr
   character*6 rank
   integer i,j,k
   integer ylen,zlen
   if(Yst==1 .and. Yed ==Ymax)then
        ylen=Yed-Yst+1
   else if( Yst==1 )then
        ylen=Yed-Yst+2
   else if(Yed ==Ymax)then
        ylen=Yed-Yst+2
   else
        ylen=Yed-Yst+3
   endif

   if(Zst==1 .and. Zed ==Zmax)then
        zlen=Zed-Zst+1
   else if(Zst==1)then
        zlen=Zed-Zst+2
   else if(Zed ==Zmax)then
        zlen=Zed-Zst+2
   else
        zlen=Zed-Zst+3
   endif

if(myid<10)then
   write (rank,150) myid
else if(myid <100)then
   write (rank,151) myid
else if(myid<1000)then
   write (rank,152) myid
else
   print*,'myid >= 1000 !'
   stop
endif
150     FORMAT (I1)
151     FORMAT (I2)
152     FORMAT (I3)
   write (rr,200) 100000+istep
200        FORMAT (I6)


   open(11,file='T'//TRIM(rr)//'_'//TRIM(rank)//'.dat')
   write(11,*) 'VARIABLES = X, Y, Z, U, V, W ,P'
   write(11,*) 'ZONE I=', Xmax, ', J=', ylen, ', K=', zlen, ', F=POINT'

   do k=Zst-1,Zed+1
   do j=Yst-1,Yed+1
if(k>0.and.k<=Zmax .and. j>0 .and. j<=Ymax)then
   do i=1,Xmax
              write(11,99) i*dx,j*dy,k*dz,u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k)
99            format(7E15.6)
   enddo
endif
   enddo
   enddo

!bak    WRITE(11,*) 'GEOMETRY M=GRID, FC=WHITE,F=POINT'
!bak    WRITE(11,*) '2'
!bak      WRITE(11,*) '4'
!bak      WRITE(11,*) x1,y1
!bak      WRITE(11,*) x1,y2
!bak      WRITE(11,*) x2,y2
!bak      WRITE(11,*) x2,y1
!bak      WRITE(11,*) '2'
!bak      WRITE(11,*) x1,y1
!bak      WRITE(11,*) x2,y1


   close(11)

end subroutine dumpp

subroutine dumpp_ave()
   use fangzhu
   implicit none

   character*6 rr
   character*6 rank
   integer i,j,k
   integer ylen,zlen
   if(Yst==1 .and. Yed ==Ymax)then
        ylen=Yed-Yst+1
   else if( Yst==1 )then
        ylen=Yed-Yst+2
   else if(Yed ==Ymax)then
        ylen=Yed-Yst+2
   else
        ylen=Yed-Yst+3
   endif

   if(Zst==1 .and. Zed ==Zmax)then
        zlen=Zed-Zst+1
   else if(Zst==1)then
        zlen=Zed-Zst+2
   else if(Zed ==Zmax)then
        zlen=Zed-Zst+2
   else
        zlen=Zed-Zst+3
   endif

if(myid<10)then
   write (rank,150) myid
else if(myid <100)then
   write (rank,151) myid
else if(myid<1000)then
   write (rank,152) myid
else
   print*,'myid >= 1000 !'
   stop
endif
150     FORMAT (I1)
151     FORMAT (I2)
152     FORMAT (I3)
   write (rr,200) 100000+istep
200        FORMAT (I6)

   open(11,file='ave'//TRIM(rr)//'_'//TRIM(rank)//'.dat')
   write(11,*) 'VARIABLES = X, Y, Z, U, V, W ,P'
   write(11,*) 'ZONE I=', Xmax, ', J=', ylen, ', K=', zlen, ', F=POINT'

   do k=Zst-1,Zed+1
   do j=Yst-1,Yed+1
if(k>0.and.k<=Zmax .and. j>0 .and. j<=Ymax)then
   do i=1,Xmax
              write(11,99) i*dx,j*dy,k*dz,u_ave(i,j,k)/(istep-steps_st),&
        &   v_ave(i,j,k)/(istep-steps_st),w_ave(i,j,k)/(istep-steps_st),p_ave(i,j,k)/(istep-steps_st)
99            format(7E15.6)
   enddo
endif
   enddo
   enddo

   close(11)

end subroutine dumpp_ave




program simple
   use fangzhu
   implicit none

   integer i,j,k
   real(8) con_err,err

   real*8 st,ed,sst,sed

   call mpi_init(ierr)
   call MPI_COMM_SIZE(MPI_COMM_WORLD,nps,ierr)
   call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)

   call init()

   call alloc_mem()

   call set_bdy()

   call set_init_flux()

       if(myid==0) print*,'end set_init_flux !'
        call flush(6)


  istep=0

        counter=0

   do istep=1,MaxStep
        con_err=0.0
        err=0.0

        call update(u)
        call update(v)
        call update(w)
        call update(p)

   ui=u
   vi=v
   wi=w
   pi=p
   uf0=uf
   vf0=vf
   wf0=wf
   pf0=pf
   mf0=mf
    


        call sub1_momentum(flagd11,flagd12,flagd13)

        call flush(6)

        call sub3_pressure(flagd11,flagd12,flagd13)

        call flush(6)

        do k=Zst,Zed
        do j=Yst,Yed
        do i=1,Xmax   
           con_err=max(con_err,abs(mf(i+1,j,k,2)-mf(i,j,k,2)+mf(i,j+1,k,1)-mf(i,j,k,1)+mf(i,j,k+1,3)-mf(i,j,k,3)))
        end do
        end do
        end do
      CALL MPI_Allreduce(con_err,err,1,MPI_DOUBLE_PRECISION,MPI_Max,MPI_COMM_WORLD,ierr)
        
        counter=counter+1
       if(myid==0) print*,'istep=',istep,'counter=',counter,'con_err=',err

     ui=u
     vi=v
     wi=w
     pi=p
     uf0=uf
     vf0=vf
     wf0=wf
     pf0=pf
     mf0=mf


    if(istep>steps_st)then
        do k=Zst-1,Zed+1
        do j=Yst-1,Yed+1
      if(k>0.and.k<=Zmax .and. j>0 .and. j<=Ymax)then
        do i=1,Xmax
           u_ave(i,j,k)=u_ave(i,j,k)+u(i,j,k)
           v_ave(i,j,k)=v_ave(i,j,k)+v(i,j,k)
           w_ave(i,j,k)=w_ave(i,j,k)+w(i,j,k)
           p_ave(i,j,k)=p_ave(i,j,k)+p(i,j,k)
        end do
      endif          
        end do
        end do

     if (mod(istep,NUMAX)==0) then 
        call dumpp_ave()
     end if
        
    endif

     if(istep<=100)then
        NUMAX=50
     else
        NUMAX=200
     endif
     if (mod(istep,NUMAX)==0) then 
        call dumpp()
     end if

   enddo


   call dealloc_mem()

   if(myid==0) print*,'success!'

   call mpi_finalize(ierr)

   stop


end program
